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ABSTRACT 

In a radiatively heated and cooled medium, the thermal instability is a plausible mechanism for 
forming clouds, while the radiation force provides a natural acceleration, especially when ions recom¬ 
bine and opacity increases. Here we extend Field’s theory to self-consistently account for a radiation 
force resulting from bound-free and bound-bound transitions in the optically thin limit. We present 
physical arguments for clouds to be significantly accelerated by a radiation force due to lines during 
a nonlinear phase of the instability. To qualitatively illustrate our main points, we perform both 
one and two-dimensional (I-D/2-D) hydrodynamical simulations that allow us to study the nonlin¬ 
ear outcome of the evolution of thermally unstable gas subjected to this radiation force. Our I-D 
simulations demonstrate that the thermal instability can produce long-lived clouds that reach a ther¬ 
mal equilibrium between radiative processes and thermal conduction, while the radiation force can 
indeed accelerate the clouds to supersonic velocities. However, our 2-D simulations reveal that a 
single cloud with a simple morphology cannot be maintained due to destructive processes, triggered 
by the Rayleigh-Taylor instability and followed by the Kelvin-Helmholtz instability. Nevertheless, the 
resulting cold gas structures are still significantly accelerated before they are ultimately dispersed. 

Subject headings: hydrodynamics — instabilities — radiative transfer — methods: numerical 


1. INTRODUCTION 

The thermal stability of astrophysical gases has been 
extensively investigated, both analytically and numeri¬ 
cally, in many different physical contexts, e.g., in star 
formation regions, in gas accreted by forming galaxies, 
and in interstellar and intergalactic media. The phys¬ 
ical basis for gas in thermal equilibrium to be subject 
to thermal instability (TI) was introduced in a classic 
paper by Field (1965), who developed the linear theory 
accounting for thermal conduction. Later Balbus (1986) 
generalized Field’s theory by relaxing the equilibrium as¬ 
sumption and permitting the background flow to be time- 
dependent. 

One specific instance where gas can be prone to Tl 
is when both cooling and heating are due to radiative 
processes. This situation occurs for example in active 
galactic nuclei (AGN) where, for a temperature of about 
10^—10® K, the cooling is dominated by free-free emission 
while the energy gain is dominated by X-ray absorption 
and Compton scattering of energetic photons produced 
by an accreting black hole (e.g., Krolik et al. 1981, Lepp 
et al. 1985). In this situation, the radiation field is al¬ 
most in the free-streaming limit and consequently the gas 
will receive a non-zero net push in a direction away from 
the source of radiation. Our objective here is to demon¬ 
strate that this push can be substantial due to a large 
increase in the bound-free and bound-bound opacities of 
the cold gas during a nonlinear phase of the TL 
In this paper, we present results from our detailed 
study on the nonlinear dynamics of thermally unstable 
gas that is accelerated by the same photons that estab¬ 
lish its thermal environment. In §2, we describe the the¬ 
oretical underpinnings on which this paper is premised. 
In §3, we outline the equations solved and the explicit 
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form of both the energy and momentum source terms. 
In §4, we present our methods and results from 1-D and 
2-D simulations. We finish in §5 where we discuss our 
findings. 

2. THEORETICAL MODEL/EXPECTATIONS 

To facilitate a simple comparison between our simula¬ 
tions and Field’s theory (i.e. the linear growth rates), 
we include thermal conduction, assume the gas to be ini¬ 
tially stationary, homogenous, and non-magnetized, and 
we do not include gravity. We model the formation and 
evolution of optically thin clouds whose thermal equilib¬ 
rium state is controlled by impinging radiation by adopt¬ 
ing a realistic prescription for heating and cooling appro¬ 
priate for gas in an AGN environment. Specifically, we 
use a net cooling function, C, that includes the following 
radiative processes: (I) Compton and inverse-Compton 
scattering, (2) photoionization and recombination, (3) 
Bremstrahlung, and (4) line-emission. In an optically 
thin case, C depends on the gas temperature, T, and 
photoionization parameter, ^ = AiriFx/n, where iPx is 
the X-ray flux and n is the gas number density}^ 

The following theoretical picture forms the basis of our 
expectations and motivates our simulation setup. For a 
given gas in equilibrium at temperature T^q will have 
heating balancing cooling [i.e., C{Teq,^) = 0]. There¬ 
fore, it will occupy one point on the radiative equilib¬ 
rium curve (the solid line which is hereafter denoted the 
S-curve) shown in the top panel of Fig. I. Suppose that 
initially this point is at a stable location on the S-curve 
such as at position 1 in Fig. 1, but some physical event 
transpires that results in a reduction of iPx so that the 
gas finds itself out of equilibrium at location 2 with, for 
example, ^2 = 190- It will quickly cool to reside on the 
S-curve at location 3 (with ^3 = ^2)- However, at this 

^ The units for ^ are erg cm s~^ and we do not cite them 
throughout the remaining part of this paper. 
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Fig. 1. — Temperature and opacity dependence on the pho- 
toionization parameter expected in an AGN environment. Top 
panel: The solid line is the S-curve found by solving £(T, = 0. 

The region above (below) this line is patterned light blue (red) 
to denote cooling (heating), as gas in this region is above (be¬ 
low) the equilibrium temperature. The dashed curve (defined by 
\dCIT')ldT\-p = 0) marks the isobaric instability criterion. Ther¬ 
mally stable gas must lie below this curve; gas anywhere in the grey 
region cannot settle on the S-curve in this region without being un¬ 
stable. The dotted line shows a constant pressure slope. Stable gas 
at location ( 1 ) is envisioned to be suddenly subjected to a reduced 
flux, placing it at location (2), where it is unstable. This gas will 
rapidly cool (nearly isochorically) until it reaches a new thermal 
equilibrium which is now unstable [marked as location (3) on the 
equilibrium curve]. We begin our simulations at this new equilib¬ 
rium to follow the growth of an isobaric perturbation. As one can 
anticipate, the perturbation grows maintaining pressure equilib¬ 
rium even during the non-linear phase and the points representing 
gas move along the dotted line; in particular those representing 
the cold gas move toward yet another thermal equilibrium location 
(4) which is now stable. All points in the computational domain 
of 1-D run RFLDX (radiation force due to X-rays and lines) are 
over-plotted as blue (T < Teq) and red (T > Teq) dots to indicate 
the final state of the gas. Bottom panel: Gas opacity in the 
units of the Thomson opacity as a function of The solid red line 
represents bound-bound opacity, Mmax- This opacity can become 
orders of magnitude larger than bound-free opacity ax (shown as 
the solid black line). The solid blue line represents the opacity of 
the most opaque line. To gauge the increase in opacity for the 
cloud formed in the top panel [location (4)], the two vertical lines 
mark the initial conditions (^3 = 190) of the gas and the eventual 
location of the cloud core (^4 p:: 73). 

position, gas is thermally unstable to perturbations with 
constant pressure p, as it violates Field’s criterion for 
stability, [dC/dT]p > 0. More generally, only the re¬ 
gion below the dashed line is thermally stable according 
to Balbus’ criterion [d{C/T)/dT]p > 0, with the con¬ 
sequence that gas occupying the shaded grey region in 
Fig. 1 will evolve to points on the S-curve outside of the 
grey region. Continuing with our example, the slight¬ 
est density perturbation present in the gas at location 3 
will grow exponentially at the linear theory rate and at 
nearly constant pressure until the perturbation becomes 


nonlinear. Rapid nonlinear growth (still at nearly con¬ 
stant pressure) will then commence (e.g., Burkert & Lin 
2000 ) and result in the formation of a cloud with a core 
density approximately determined by the intersection of 
the dotted line and the S-curve (position 4), which is 
where the gas in the core is stable. 

The bottom panel of Fig. 1 illustrates how various gas 
opacities depend on Solid black and red lines show the 
bound-free opacity, cr^, and the total line opacity, Mmax, 
respectively. The solid blue line represents the opacity 
of the most opaque line. One can see that for ^ 3 , the 
bound-free and bound-bound opacity is negligible. (The 
vertical line in the grey region marks the location ^3 = 
190.) However, once the thermally unstable gas forms a 
dense cloud, ax and Mmax are significantly increased, as 
indicated by the left vertical line at ^4 « 73. Therefore 
the cloud can be accelerated by the same source that 
heats it. 

3. GOVERNING EQUATIONS 

To confirm and quantify our expectations, we numeri¬ 
cally solve the equations of hydrodynamics: 

-f V • (pv) = 0, (1) 

+V ■ {pVV +p\) =frad, ( 2 ) 

f) F 

— + V • [(£; + p)v] = -pC + «:o V^T. (3) 

In the above, p is the gas density, v is the fluid velocity, 
p is the gas pressure, I is the unit tensor, frad is the 
radiation force (see below), and E = e + pv^ 12 is the 
total energy density, with e the internal energy density. 
This system is closed with an ideal gas equation of state, 
e = pfipi — l)j using an adiabatic index 7 = 5/3. We 
estimate the conduction coefficient kq using the value 
for a fully ionized plasma (Spitzer 1962) evaluated at 
the equilibrium temperature and density of position 3 
in Fig. 1 (hereafter a suffix ‘eq’ denotes evaluation of a 
quantity at [^ 3 , Teq]). 

We denote the functional dependence of the net cooling 
function as £ = A — F, with 

A = + Lbb) [ergg“^ s“^], and (4) 

firup 

r=— (Gc + G^) [ergg-^s-^], (5) 

finip 

where p is the mean molecular weight (set to 1.0 in this 
work) such that n = prupp (with nrip the proton mass), 
Lff and Lbb are the cooling rates due to free-free and 
bound-bound transitions, and Gc and Gx are the heat¬ 
ing rates due to Compton and X-ray heating, respec¬ 
tively (all four rates are in units of erg cm^ s“^). For 
Lff and Gc^ we use the well-known analytic formulas 
based on atomic physics; see the appendix. For Lbb 
and Gjjf) meanwhile, we use the analytical fits given by 
Blondin (1994), who found good agreement (to within 
25%) between his approximate rates and those result¬ 
ing from detailed photoionization calculations. Blondin’s 
formulae, also provided in the appendix, assume an op¬ 
tically thin gas of cosmic abundances illuminated by 
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a Tx = lOkeV/ke bremsstrahlung spectrum (with fcs 
Boltzmann’s constant). 

To evaluate frad, we consider the source of radiation to 
be located at a: = — oo, giving 

n P^tot^tot ^ r —2 —2l ir‘\ 

W =-a: gem s , ( 6 ) 

c 

where Ttot is the total continuum flux, c is the speed 
of light, and a tot is an effective total opacity, individual 
contributions to which can be self-consistently estimated 
as we describe below. 

The X-ray flux is set by the photoionization param¬ 
eter Therefore, we can specify the product Otot^tot 
by assuming that there is nonzero continuum opacity 
only for X-rays and nonzero line opacity for UV pho¬ 
tons. With this assumption we parametrize the UV flux 
J^uv in terms of the X-ray flux using /yy = •^uv/-^x 
and use the following expression to estimate the opacity 
and flux product: 


^tot^tot — 


(1 + /uv) + Crx + f\JvMyna.x ^X- (7) 


Here, Ue is the mass scattering coefficient of the free 
electrons. The term (1 -|- /uv) in equation Q is due 
to electron scattering, while the term ux is an effec¬ 
tive X-ray opacity in units of We compute ax as 
(dTrGx./i/O/CTe, where Gx,h is the heating part of Gx 
(see the appendix). The last term in equation 0 , Vf^iaxj 
is the maximum force multiplier characterizing the line 
contribution to the total opacity. We evaluate M^ax fol¬ 
lowing Stevens & Kallman (1990, SK90 hereafter) with 
a modification, as described below. Note that in AGN, 
/uv > 1 - 

The maximum force multiplier parametrizes the in¬ 
crease in the scattering coefficient due to spectral lines 
when all of the lines are optically thin. In this opti¬ 
cally thin limit, the multiplier is just a sum of opac¬ 
ity contributions from all the lines and depends only on 
the gas composition, ionization, excitation and oscillator 
strengths [e.g., see eqs. 10 and 11 in Castor, Abbott, 
& Klein, (1975); CAK hereafter]. Following Owocki, 
Castor, & Rybicki (1988), who used a modified CAK 
method, one can evaluate this maximum multiplier as 

Afmax = AcAk(1 ^)Vniaxi (^) 

where fccAK is proportional to the total number of lines, 
a is the ratio of strong to weak lines, and finally r?max = 
^L.max/o’e is a dimensionless measure of opacity of the 
most opaque line (with KL,niax being the line opacity co¬ 
efficient of the thickest line). 

SK90 carried out detailed photoionization calculations 
for a radiative environment appropriate for X-ray sources 
and parametrized their results in terms of the above ex¬ 
pression for Mmax by allowing kcAK to be T-dependent 
and ?7max to be /-dependent. Instead of the fit for 
fccAK('T) from SK90, we use equation (17) of Proga 
(2007) due to Kallman (2006, private communication), 
as this may better represent the increase in the number 
of lines with decreasing temperature in AGN. This ex¬ 
pression and that for ?7max(/)j which is equation (19) of 
SK90, are provided in the appendix. Both of these fits 
were generated assuming a = 0.6. In the bottom panel 
of Fig. I, we plot ax along with Mmax- Notice that 
can be roughly a few thousand for gas ionized by 


a weak radiation field, whereas it decreases asymptoti¬ 
cally to zero for highly ionized gas as the radiation field 
becomes stronger. 

4. METHODS AND RESULTS 

We solve equations 0-0 in 1-D and 2-D using the 
CTU integrator. Roe flux, and explicit conduction mod¬ 
ule of the MHD code Athena (Stone et al. 2008). We 
modify the original version of the code by adding the 
momentum and heating and cooling source terms in the 
same way that Athena’s built-in static gravitational po¬ 
tential source term is implemented to achieve 2 nd order 
accuracy in both space and time. We use a less ac¬ 
curate method for integrating the conduction term in 
time in our 2-D simulations than in our 1-D simula¬ 
tions. Specifically, we use Athena’s super time-stepping 
scheme (STS; see O’Sullivan & Downes 2006), although 
we note that a 2nd order accurate in time STS algorithm 
does exist (Meyer et al. 2012) and we are testing it for 
future use. 

4.1. Initial and Boundary Conditions 

Given the atomic physics behind our S-curve and the 
opacities in Fig. 1, the following free parameters govern 
our problem: the wavenumber and density amplitude of 
the TI perturbation, k and Sp, as well as the ratio of the 
initial sound crossing and thermal times Uh/tsc, which 
together determine the number of clouds and their for¬ 
mation time; the initial photoionization parameter /a, 
which controls the intensity of the radiation field and 
the equilibrium temperature Tg^; the equilibrium pres¬ 
sure, namely the product UegT^q, which sets the physical 
units of the cloud and its environment; and Anally fuv^ 
which parametrizes the shape of the spectral energy dis¬ 
tribution in a simple way. 

Our initial conditions are full wavelength profiles of the 
TI condensation mode, found from equations (11)-(14) 
in Field (1965), applied to density, velocity, and pres¬ 
sure. We adopt n^gTeq = 10^^ K cm“^ in accordance 
with AGN observations and their modeling (e.g., David¬ 
son & Netzer 1979; Krolik et al. 1981). The length scale 
of the perturbation is fixed by the adiabatic sound speed 
at position 3 in Fig. 1, Cgq, and a choice for the initial 
sound crossing time, Gc- We chose Gc equal to the ini¬ 
tial thermal time, tth = i^eq !which results in 
near maximum linear growth rates for the condensation 
mode. These rates are obtained by solving the dispersion 
relation in Field (1965; eq. 15). For both 1-D and 2-D 
simulations, we use periodic boundary conditions and set 
the domain size in the x-direction, Ax, equal to the per¬ 
turbation wavelength = Ceqtsc- The amplitude of the 
density perturbation is Sp = 5.0 x p^q. 

With this setup, the unstable region of Fig. 1 is 
parametrized entirely by /a, and varying this parame¬ 
ter leads to the formation of clouds with substantially 
different properties. A realistic value for AGN is likely 
/a = 500, as this results in a pressure photoionization 
parameter S = {IFx/c)/{ueqkBTeq) ~ 9.0, and the AGN 
environment is expected to be hospitable to clouds for 
S < 10 (e.g., Krolik et al. 1981). For /a = 500, the 
linear growth rate of the TI is comparatively small, tak¬ 
ing more than 400 days for the density of the cold gas 
to double, and estimates based on Fig. 1 indicate that 
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a 1-D (or 2-D planar) cloud will form at ~ 20 with 
a width If. ^ 0.5 AU, a temperature ~ 4 x 10^ K, a 
number density ric 3 x 10 ®cm“^, and a density con¬ 
trast of X ~ 30. The radiation force due to lines can be 
very powerful, with M^ax at least 10 ^. 

This more realistic cloud is, however, very optically 
thick for many strong UV lines. (Moreover, as we dis¬ 
cuss in §5, it is challenging to accurately resolve in multi¬ 
dimensions.) Our present simulations are designed to ex¬ 
plore the optically thin regime, as this regime allows us 
to focus on the purely hydrodynamical effects of cloud 
formation and acceleration without the further compli¬ 
cations involved when solving the equations of radiation 
hydrodynamics (RHD; cf. Proga et al. 2014). This re¬ 
stricts ^3 to a very narrow range of values corresponding 
to larger, less realistic values of S, as there is obviously an 
upper limit on the density of clouds whose acceleration 
we can accurately model without RHD. 

To estimate this upper limit on the density, we as¬ 
sume the cloud forms with both a constant density 
and width (which will be seen to be very nearly the 
case here), so that we can write Tymax = T^.-max/Ts, where 
TL.max = K.L.maxpJ'c IS Optical depth of the thickest line 
and Ts = (TepJc is the electron scattering optical depth of 
the cloud. Demanding that the cloud to be optically thin 
to all bound-bound transitions (i.e. Tj^^max < 1 ) requires 
VmaxTs < 1, Or 

Pc ^ i^e^cVmax) • (9) 

We chose the value ^3 = 190 since it produces a cloud 
with the highest density contrast that satisfies this in¬ 
equality at all times in 1-D. (In §4.4, we verify our opti¬ 
cally thin assumption in both 1-D and 2-D using a more 
accurate estimate.) For ^3 = 190, Mmax does not exceed 
40. To explore the effects of the stronger line force, we 
set f^v = 10 , which is guided by observational results 
from Zheng et al. (1997) and Laor et al. (1997). 

4.2. Simulations 

We have performed over 30 simulations exploring a 
variety of parameters and numerical setups. Here we 
report in some detail on a set of four simulations that 
differ only by the applied radiation force: RF (electron 
scattering only), RFX (electron scattering plus X-ray ab¬ 
sorption), RFLD (electron scattering plus line-driving), 
and RFLDX (electron scattering, line-driving, and X-ray 
absorption). We will especially focus on run RFLDX, as 
this case is most representative of the physical condi¬ 
tions in AGN. Compared to a more realistic AGN cloud 
described above, for our adopted value of ^3 = 190 the 
cloud growth is much faster, taking only 2.4 days to 
double in density, while 4 ~ 6.5 x 10 ^° cm (0.004 AU), 
Tc « 7.0 X lO^K, ric « 1.4 x 10® cm-^, y « 8.0, and 
S « 19. The physical units corresponding to our numer¬ 
ical results are listed in Table 1. 

4.3. Results of 1-D Simulations 

Despite the different radiation forces, the clouds in all 
four runs are formed at the same time and with the same 
density and temperature contrasts, and the gas therefore 
traces the same ‘tracks’ on the T — f plot in Fig. 1. The 
over-plotted red and blue dots in Fig. 1 show the tracks 
for RFLDX at t = 120 tsc, which represents the time 


TABLE 1 
Physical units 

(corresponding to ^ = 190 & tth/tsc = 1) 


Quantity 

Value (cgs units) 


peq 

8.65 

X 

10" 

17 

' g cm 

-3 


5.17 

X 

10" 

cm~^ 


T 

eq 

1.93 

X 

10^ 

K 


Ceq 

5.16 

X 

10® 

cm s~^ 


tsc 

6.03 

X 

10"^ 

s 



3.11 

X 

lOi 

^ cm 


Tx 

7.82 

X 

10® 

erg 

cm-2 

m 

1.58 

X 

10" 

erg 


Agg 

3.97 

X 

10® 

erg g-" 

s-l 


where the flow reached a thermal steady state (see below 
for more details). Note that the red tracks do not reach 
the radiative equilibrium curve (i.e. £ = 0 ), but rather 
an equilibrium curve given by pC = nfS/^T. 

Also note that there are tracks occupying an unstable 
(according to Balbus’ criterion) region in Fig. 1, namely, 
the tracks within the grey region that are above the 
dashed line. Given that the gas in the cloud core oc¬ 
cupies location 4 and is in pressure equilibrium with the 
medium, it must be the case that some portion of the 
gas occupies this unstable region in order for the density 
and temperature to be continuous everywhere. These 
‘unstable’ tracks correspond to the gas in the conductive 
interfaces of the cloud. In Fig. 2, we plot profiles of the 
solution overplotted in Fig. 1, and the width of the in¬ 
terfaces can be judged from the density panel. The local 
Field length in the interfaces is close to the initial equi¬ 
librium value, \f = 2'KyJ KeqTeq! {peq^-eq) Ax « 0.19 A^,, 
while the interface width is two or three times smaller 
than this. Gas is permitted to be thermally unstable at 
regions smaller than the Field length when the stabiliz¬ 
ing influence of conduction (measured by the heat flux 
that is shown in the fifth panel of Fig. 2) is large enough. 

We hnd that in all of these runs the gas arrives at 
a state of thermal equilibrium by t = 90tsc- The sec¬ 
ond from bottom panel in Fig. 2 shows how this equilib¬ 
rium state is possible. Both the net cooling function (red 
curve) and the conduction term (black curve) are posi¬ 
tive at the interfaces of the cloud and negative in the hot 
medium. These terms are of opposite sign in equation ([^ 
and therefore balance each other. The compression term 
(cyan) is negligible at this time, but it was the dominant 
term when the cloud was forming. 

The radiation force prevents the gas from reaching a 
mechanical equilibrium state. Rather, in each case the 
cloud core undergoes dynamical changes (i.e. the pres¬ 
sure and velocity profiles adjust) to permit nearly uni¬ 
form acceleration. To show this, we plot the net flow 
acceleration in the bottom panel of Fig. 2 (cyan line), 
which is the sum of the other curves displayed. Line driv¬ 
ing operates almost uniformly throughout the cloud core. 
As can be seen by either the acceleration or the pressure 
panel, the response of the gas pressure is to exert a nearly 
constant drag force on the cloud to compensate for the 
driving force, while the medium is pushed along since it 
has nowhere else it can go in 1-D. The adjustment of the 
forces is obtained shortly after the cloud is formed, with 
the acceleration profiles resembling those shown in the 
bottom panel at around t = 55 tsc for runs RFLD and 
RFLDX. Some prohles (especially velocity) continue to 
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Fig. 2.— Profiles of run RFLDX in 1-D at time 120£sc- The resolution is Nx = 1024 zones. Numbers above panels are offsets, e.g., the 
velocity ranges from about 1.3843 Ceg to 1.3851 Ceg. The dotted vertical line indicates the position of the maximum density gradient. The 
2nd from bottom panel compares the heating and cooling rates R in the energy equation, while the bottom panel compares the various 
accelerations. Solid (dashed) portions of lines in these panels indicate positive (negative) values, e.g. conduction transfers heat into the 
interfaces at the expense of the medium, while the specific pressure force points to the left (opposite the cloud motion) in the cloud core 
and to the right elsewhere. 


undergo changes until t = 120 tsc; the shapes shown in 
Fig. 2 are maintained as the cloud continues to accelerate 
to high Mach numbers. 

The results from our 1-D simulations confirm our ba¬ 
sic picture for cloud formation and acceleration. As a 
next step, we perform 2-D simulations where destructive 
processes may change our results. Our primary concern 
is the Rayleigh-Taylor (RT) instability, and subsequently 
the Kelvin-Helmholtz (KH) instability. Fig. 2 shows that 
the left interface is RT stable, as the density increases in 
the direction of acceleration. However, the right interface 
is likely unstable due to the adverse density arrangement 
(e.g., Krolik 1977, 1979; Mathews & Blumenthal 1977; 


Jacquet & Krumholz 2011; Jiang et al. 2013). 

4.4. Results of 2-D Simulations 

The simplest extension of our 1-D simulations to 2-D 
is to form a cloud with a planar slab configuration. The 
initial conditions and overall setup is as before. However, 
to fully explore 2-D effects, we now break the uniformity 
in the y-direction by introducing a perturbation with a 
wavelength the size of the domain in the y-direction (i.e. 
Aj, = Aa;/2) and Sp = 5 x p^q. 

In Fig. 3, we present a comparison of our four 1-D 
runs and their 2-D counterparts, and we also verify our 
optically thin assumption. The maximum density of the 
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Fig. 3. — Comparison of all runs in 1-D and 2-D. Dashed (solid) 
lines denote 1-D (2-D) runs. In both panels, all curves nearly over¬ 
lap during the nonlinear cloud formation process, which completes 
at time 5^ 50 tsc- In the top panel, we also verify that the cloud 
in run RFLDX is optically thin to UV radiati on b y calculating an 
estimate to the optical depth using equation (|10[ ). This estimate 
in 1-D (2-D) is given by the dashed (solid) black line. The dot¬ 
ted vertical lines mark the times corresponding to the snapshots in 
Fig. 4. 

cloud versus time is plotted in the top panel. It is clear 
that there is no significant difference in any of the runs 
during the cloud formation process, which ends at t « 
50tsc- As mentioned in §4.1, for runs RFLD and RFLDX 
to be optically thin to UV photons, we require t/, ^lax < 
1. We estimate as|^ 

TL.max = CTe / p{x)dx, (10) 

Jo 

where Tymax.go denotes only those values of ?7max in the 
range [0.9 r^max, ?7max]- We use this range to be able 
to identify gas at a constant opacity in our numerical 
representation of rymax- The black curves in the top 
panel of Fig. 3 show this estimate for max in both 1-D 
(dashed line) and 2-D (solid line). In the latter calcula¬ 
tion, we consider a ray through the center of the cloud at 
y = 0.25Aa;. Overall, the cloud can indeed be considered 
optically thin at all times, except possibly at its center 
during a very short period of the acceleration phase in 
2-D, which as will be made clear below, coincides with 
when the cloud is significantly lengthened by the onset 
of disruptive processes. 

The bottom panel of Fig. 3 shows the average velocity 
of the cloud versus time, where we define the cloud as 

^ Note that the ‘expanding’ optical depth formula, r/^^max = 
fTe r;max (see CAK) is not valid here because the 

Sobolev length vth/\d‘V/dx\ is much greater than the density scale 
height. That is, the scale over which the velocity changes by the 
thermal width of the line ~ 20 kms“^) is much greater than 
the interface width (~ Af), the scale over which the density (and 
hence opacity) changes appreciably. 


being the gas to the left of the grey region in Fig. 1 (i.e. 
^ < 121.2, which corresponds to p > 1.57peq)- It indi¬ 
cates that significant acceleration only takes place after 
the cloud formation process has ended, as line opacity 
is only activated once the cloud forms. While 1-D runs 
of RFLD and RFLDX uniformly accelerate to supersonic 
speeds, the acceleration is suddenly halted in 2-D around 
t « 66 tsc- 

In Fig. 4, we plot snapshots of run RLFDX, to illus¬ 
trate the very different fate of a rapidly accelerated cloud 
in 2-D. The first frame shows that the slab formed at 
t = 45tsc- The initial perturbation in the y-direction 
grows into a slight over-density region in the center of the 
slab, which then undergoes greater acceleration than its 
surroundings and causes a small bulge to appear around 
t « 55tsc (not shown). Viewing this bulge as a pertur¬ 
bation along the surface of the slab, the basic criteria 
for the RT instability is satisfied: heavy fluid is pushing 
against light fluid. The same conclusion applies to runs 
RFX and RFLD, although for run RFX it will take much 
longer (several hundred tsc) for the bulge to grow due to 
the weak acceleration. Run RF, however, which has a 
constant acceleration due to Thomson opacity, evolves 
identically in 1-D and 2-D for all time; any density per¬ 
turbation in the y-direction receives the same push as 
any other point in the flow. 

The remaining frames in Fig. 4 reveal how the breakup 
of the cloud ensues as the RT instability develops and 
soon becomes accompanied by the KH instability. First 
the bulge evolves to become mushroom-shaped, forming 
a structure resembling the classic RT ‘spike’, which fea¬ 
tures prominent KH ‘rolls’ at t = 66tsc made possible 
by the increased relative velocity between the cloud and 
medium. The halting of the acceleration happens around 
this time. The connecting plume then disperses (i.e. its 
gas is heated) as the spike separates further from the slab. 
The slab would likely disperse also, due to the mass lost 
to the spike, but instead it is somewhat thickened by the 
approach of the spike from the backside. This collision 
perspective is shown in the final panel at t = 75tsc- 

The cloud is eventually either completely dispersed 
back into the medium, or coexists with it in a disordered 
manner in what could be called a clumpy flow once the 
vertical symmetry is lost. The simulations presented here 
do not let us make any definitive statements because we 
noticed that our solutions become resolution dependent 
at about time 73 tsc for run RFLDX. This loss of conver¬ 
gence is shown in Fig. 5, where for four different resolu¬ 
tions, we plot the mass fraction of the three components 
of the gas: (i) the cloud (blue), again defined as gas with 
^ < 121.2 or p > 1.57 peq (ii) the interface (black), de¬ 
fined as the portion of the gas in the grey region that is 
unstable, i.e. the tracks above the dashed line in Fig. 1 
with 121 < ^ < 278 or 0.68 Peo < P < 1-57 Peq] (in) the 
medium (red), defined as the (stable) gas with ^ > 278 
or p < 0.68 Peq- The mass fraction is defined as the mass 
of each component divided by the total mass and is given 
by m = {NxNy)~^Y,ij Pij!Peq, where and Ny are the 
number of grid zones in the x and y directions and the 
sum ranges over all zones (i,j) that satisfy one of the 
criteria (i)-(iii). Once the mass fractions for resolutions 
Nx = 1024 and Nx = 2048 differ, we cannot claim to 
accurately follow the cloud’s evolution. 
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0.67 


0.36 


2.4 


Fig. 4. — Density snapshots of run RFLDX in 2-D in units of p^q. The domain size is [A^, Ai,/2] with resolution A^y] = [1024,512]. 
Since the cloud continually advects through the domain boundaries, the images are manually aligned for visual comparison. Velocity arrows 
are overlaid after subtracting the mean x-velocity of the cloud (displayed in the upper right corner) from v^, the cloud being defined as gas 
with p > l.STpeg. Time is shown in the lower left corner of each panel. 


The mass fractions provide a complementary descrip¬ 
tion of the cloud evolution depicted in Fig. 4. The cloud 
appears fully formed by t « 50tsc, which coincides with 
the peak mass fraction of the medium, but mass keeps 
piling on until t « 63tsc- Indeed, we observe that the 
velocity arrows at t = 50 <sc in Fig- 4 point toward the 
cloud, indicating continued growth at the expense of the 
medium. The overall fraction of gas occupying interface 
regions is a minimum at t « 60 tsc, and this is despite the 
overall increase in the size of interface region (due to the 
bulge) because the interfaces are narrower. The cloud 
mass fraction reaches a maximum at t « 63tsc when 
the RT spike has become fully nonlinear, and thereafter 
monotonically decreases, with the mass being taken up 
entirely in the interfaces, until the RT spike becomes 
detached from the slab. During this time, the medium 
continues losing mass to the interfaces. The loss of con¬ 
vergence is likely due to the appearance of small scale 
structures as the cloud is disrupted. 

4.5. 2-D Convergence Study 

The resolution dependent late-time dynamics encoun¬ 
tered above warrants demonstrating that our results are 
converged at earlier times, especially since an accelera¬ 


tion scheme such as STS was found to be necessary in 
2-D due to the strict time constraint imposed by the CFL 
condition (At cx Aa:^ due to conduction) combined with 
the long duration of the runs. To show that our 2-D re¬ 
sults using STS are converged, in the top panel of Fig. 6 
we plot Li errors for run RFLDXa at t = 45 tsc- This run 
is a modified version of run RFLDX where we suppressed 
the contribution of the force due to electron scattering by 
a factor of 11 (by setting the term (1 -|- fuv) in equation 
0 equal to 1) to minimize the effects of advection errors 
not associated with the dominant term fuvMmax- 
For a uniform mesh the Li error norm is the quantity 
iNa;Ny)-^J2ij where Rij = \qij - qref \ is the local 
absolute error in some quantity q at the location of zone 
(i,j) with respect to some reference solution qref evalu¬ 
ated at the same location, and the sum ranges over all 
zones in the domain. As a reference solution we used the 
highest resolution simulation we could reasonably afford, 
namely [Nx,Ny] = [2048,1024]. The dashed line in Fig. 2 
marks the slope corresponding to a convergence rate of 
2 , showing that the solutions for density and velocity 
converge to the reference solution in line with Athena’s 
2nd order accurate spatial discretization. However, pres¬ 
sure shows a reduced convergence rate, which may be 







































































Fig. 5. — Mass fractions for run RFLDX in 2-D at four different 
resolutions. Blue, black, and red curves track the fraction of the 
gas contained in the cloud, the interfaces, and the medium, respec¬ 
tively; our method of differentiating these regions is described in 
§4.4. The dotted vertical lines mark the times corresponding to 
the snapshots in Fig. 4. The mass fraction of the cloud monoton- 
ically increases until the RT spike becomes fully nonlinear, after 
which it begins to monotonically decrease until the spike becomes 
a detached structure. At this point our results become resolution 
dependent. 



Fig. 6 .— Convergence study of run RFLDXa in 2-D, showing 
Li errors for resolutions 256, 512, and 1,024 in the x-direction, 
with respect to the Nx = 2,048 reference solution at time 45 f sc- 
The dashed-dotted and dashed lines show slopes corresponding to 
1st and 2nd order convergence, respectively. Pressure only has a 
convergence rate of 1, which is likely due to the reduced temporal 
accuracy of the STS scheme used to integrate the conduction term. 


due to the fact that the STS scheme is only 1st order ac¬ 
curate in time. Interestingly, we observed that reducing 
the Courant number by a factor of 6 will result in pres¬ 
sure also showing 2nd order convergence, but our results 
presented here used the default Courant number 0.8. 

In Fig. 7 we plot relative errors for pressure, defined 
as Rij/\pref\, for a horizontal slice through the domain 
at j = 0. This plot also demonstrates self-convergence 
and, not surprisingly, it shows that the largest errors are 
found at the conductive interfaces of the cloud. 

5. DISCUSSION 



Fig. 7.— Relative errors in pressure for run RFLDXa in 2-D 
at time 45 tsc- The ‘humps’ centered around x = 0.5 Aa; mark 
the conductive interfaces of the cloud, which are the most difficult 
regions to resolve and thus have the largest relative errors. 

As the next planned phase of our ongoing investigation 
of cloud acceleration initiated in Proga et al. (2014), this 
paper considered the cloud formation and acceleration 
processes simultaneously for the first time. In particular, 
we have extended the basic theory of the nonlinear out¬ 
come of TI by self-consistently solving for the dynamics 
of optically thin gas in the presence of a strong radiation 
field. Our resulting simulations have made it possible to 
study in detail the evolution of the isobaric condensation 
mode in thermally unstable gas from an initial pertur¬ 
bation to a dense, high velocity cloud. The cloud forms 
in a radiation pressure dominated environment, but the 
radiation force has practically no effect on the cloud for¬ 
mation process in either 1-D or 2-D. The reason is simply 
because there is no momentum transfer from the radia¬ 
tion to the gas unless there is also sufficient opacity, and 
the sources of opacity are not activated until the cloud 
is formed. 

The initial motivation for this work was simply to 
demonstrate that the nonlinear phase of TI leads to 
a natural mechanism to produce fast clouds via accel¬ 
eration due to lines. In so doing we were led to the 
inescapable conclusion that accelerated clouds undergo 
rapid deformation and are ultimately destroyed, thereby 
confirming long-standing assertions about the inevitabil¬ 
ity of cloud destruction (e.g., Mathews 1986; Krolik 1999 
and references therein). However, we find that optically 
thin clouds can survive long enough to be accelerated to 
relatively high velocities and travel a significant distance 
of many cloud sizes, in contrast to investigations explor¬ 
ing pre-existing optically thick clouds (Proga et al. 2014 
and references therein.) Consequently, the best hope for 
cloud-based models of AGN is to identify robust mecha¬ 
nisms for continually producing clouds (e.g., in “outflows 
from inflows” as illustrated in simulations presented by 
Kurosawa & Proga 2009 or Moscibrodzka & Proga 2013). 
It is therefore important to thoroughly study how clouds 
form via TI and how they evolve before disruptive pro¬ 
cesses take hold. 

We found that a rich set of dynamics unfolds during the 
cloud formation process. For example, for run RFLDX 
(radiation force due to X-rays and lines) in 1-D, there are 
three nonlinear phases of the TI: (i) initial growth and 
saturation (t « 40 — 46Gc); (h) evolution toward a uni¬ 
formly accelerating solution (t « 46 — 55 tgc) (hi) evolu¬ 
tion toward a thermal equilibrium state {t « 46 —90tsc)- 
Figure 4 shows that phases (i) and (ii) have both com- 
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pleted before the RT instability can develop, so these 
phases also take place in 2-D, while phase (hi) has a 
different outcome in 2-D and the uniform acceleration 
cannot be maintained. These phases will be described in 
more detail in Waters & Proga (2015, in preparation). 

Our follow up paper will also further investigate the 
growth of the RT and KH instabilities. This is an im¬ 
portant matter since the appearance of these instabilities 
governs the end phase of TI and therefore dictates the 
ultimate fate of the cloud. These instabilities develop af¬ 
ter the TI saturates. Therefore, it should be possible to 
confirm the theoretical linear growth rates by conducting 
a careful numerical perturbation analysis (e.g., by intro¬ 
ducing sinusoidal perturbations along the interfaces of a 
2-D planar slab initialized using a 1-D solution). That 
said, we did not find it possible to make a meaningful 
comparison of the growth rate of the bulge and the clas¬ 
sical RT rate in the present setup. Recalling Fig. 3, the 
bulge is already borderline nonlinear hy t = GOtgc, im¬ 
plying that the slab is still evolving thermally [i.e. in 
phase (iii) of nonlinear TI evolution] during the linear 
RT growth regime. This dynamical complication com¬ 
bined with the simplifying assumptions inherent in the 
linear theory for RT (such as constant acceleration ev¬ 
erywhere) warrant using a more controlled approach for 
isolating the development of the individual instabilities. 

The numerical setup presented here can be used for 
several different exploratory studies of cloud formation 
and acceleration. We chose the initial perturbation am¬ 
plitude Sp in the y-direction to be 10“^ that in the x- 
direction in order to arrive at a slab configuration. With 
an equal ratio (and equal growth rates), a round cloud 
will be formed in 2-D instead. Multiple clouds can be 
formed using higher wavenumber perturbations. Cloud 
fragmentation can be studied by decreasing the ratio 
tth/tsc- Classical evaporation (Cowie & McKee 1977; 
Balbus 1985) can be explored by arranging for the Field 
length to exceed the size of the cloud. This large range 
of initial configurations would be difficult or impossible 
to construct otherwise, which illustrates an obvious ad¬ 
vantage of studying cloud acceleration via the formation 
process, namely that the internal gas dynamics is self- 
consistently treated. Indeed, we find that pressure equi¬ 
librium with the surrounding medium is naturally main¬ 
tained, cloud interfaces form with a width determined by 
the conductivity, the radiation and drag forces reach a 
balance so that hydrostatic equilibrium is established in 
the reference frame of the cloud, and the cloud reaches a 
thermal equilibrium state in which heating by conduction 
is balanced by line cooling, in agreement with Begelman 
& McKee (1990). Moreover, the equilibrium location on 
the S-curve largely determines the cloud density, temper¬ 
ature, and opacity before it is accelerated, while plotting 
the evolutionary tracks on the T — ^ plane has shown 
itself to be a useful tool both for understanding the time 
evolution of the cloud and for characterizing the compo¬ 
nents of the gasj^ 

Taking into consideration the numerical requirements 
involved in this study, multi-dimensional simulations will 
likely be constrained to only explore values of ^3 < 300 
for the S-curve used here, thereby limiting the overall 

^ Simulations demonstrating this can viewed online at 
WWW . physics .unlv.edu/astro/pw 15sims. html 


density contrast of clouds to y « 18. This limitation 
arises due to the need to resolve interfaces, as simula¬ 
tions of clouds formed via TI will not be converged unless 
the conductive interfaces are resolved (Koyama & Inut- 
suka 2004). Previous numerical studies using pre-existing 
clouds without thermal conduction and with unresolved 
interfaces were able to explore much higher density con¬ 
trasts such as X = 50 (McCourt et al. 2014) and even 
X ~ 10'^ (Krause et al. 2012). We found that a realistic 
K oc T5/2 would require upwards of = 4, 096 zones 
(i.e. at least three levels of refinement if adaptive mesh 
refinement is employed) even for our modest density con¬ 
trast of X ~ 8 , as would values of x ^ 30 with continued 
use of a constant conductivity. This rapid steepening of 
the interfaces with either increased x or realistic k(T) re¬ 
sults in ever smaller transition regions between the cloud 
and the medium but does not imply that the role of ther¬ 
mal conductivity becomes less important. As a conse¬ 
quence, the conduction time step would be so small at 
these resolutions that even STS schemes would become 
impractically slow, necessitating the use of implicit tech¬ 
niques. Such extensions to this work may be needed to 
assess, for example, if the the timescale for cloud destruc¬ 
tion is sensitive to the slope of the interfaces, i.e. if it 
decreases with steeper density gradients, as was found to 
be the case when a cloud is disrupted by the passage of 
a shock (Nakamura et al. 2006). 

Future work is also needed to address important as¬ 
pects of cloud formation and acceleration neglected here. 
Magnetic fields can play a dominant role in both, as they 
can prevent the gas from condensating in the first place 
(e.g., Mathews & Doane 1990), they may help accelerate 
clouds through confinement (e.g., Arav & Li 1994), and 
they can significantly effect cloud dynamics due to the 
effects of anisotropic conduction (Choi & Stone 2011). 
Other environmental factors such as the presence of Cori¬ 
olis forces, gravity, winds, shocks, etc. could be incorpo¬ 
rated to build a more global model of cloud dynamics 
for comparison with observations. However, it is cru¬ 
cial to model 3-D effects, as both compression and ex¬ 
pansion of clouds can be enhanced by a large factor, 
and this can affect the timescale for cloud destruction. 
A direct extension of this study designed to form ini¬ 
tially spherical clouds in 3-D would violate our optically 
thin assumption and therefore requires a full RHD treat¬ 
ment using methods that can accurately capture shad¬ 
ows (e.g., Davis et al. 2012; Jiang et al. 2012). Such 
methods, which have already been employed by Proga 
et al. (2014), would also allow investigation of the so- 
called line-deshadowing instability (Owocki & Rybicki 
1984) that could lead to changes in the structure and 
time-dependence of the clouds. 
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APPENDIX 

The net cooling function £ that defines our S-curve is comprised of four heating and cooling rates. The analytic 
expressions for Lq and Gc are 


Le = 


2®7re® / 2'KkBT 


Gc = 


3hmeC^ 


3mp 


Z'^qb = 3.3 X 10 [ergcm^s ^], and 


47rmpC^ 


T 




^Tx ( 1 - 4— ) = 8.9 X ( 1 - 4— ) [erg cm^s"'], 


T 


X 


( 1 ) 

( 2 ) 


where rUe and e are the mass and charge of an electron, h is Planck’s constant, Z is the ion atomic number, gs is an 
averaged Gaunt factor, and Tx = lOkeV/ke. The analytical fits from Blondin (1994) in our notation read 


Lhh = S 


1.7 X 10"^® exp I - 


1.3 X 10® 




-24 


[erg cm®s [, and 


Gx = 1-5 X 10-21^1/4^-1/2 J^i _ cm®s-®]. 


(3) 

(4) 


Here, d is a parameter introduced by Blondin (1994); setting d < 1 mimics reducing the strength of line cooling 
when relaxing his assumption of optically thin gas. We keep <5 = 1 since we assume optically thin clouds. In ^ we 
refer to the heating part of Gx, which is Gx,h = 1-5 x [erg cm®s“^]. Finally, it is important to note 

that Blondin’s photoionization calculations were revisited and independently verified using XSTAR by Dorodnitsyn 
et al. (2008) for an incident AGN power law spectrum. Their analytical fits differ from the above only by a minor 
modification to equation ([^, which Dorodnitsyn et al. (2008) report had no significant dynamical effects on their 
simulation results. 

For line-driving, we use equation (17) from Proga (2007) 


log fccAK 


-0.383 

-0.6301ogT-P 2.138 
-3.8701ogT-P 17.528 


for log T < 4 
for 4<logr<4.75 
for logr>4.75 


and equation (19) from Stevens & Kallman (1990): 

6.9 exp(0.16^°-^) 


for 


log Tjraax — 


9.1 exp(—7.96 X 10 4^) for 


log^ < 0.5 
log^ > 0.5. 


(5) 


( 6 ) 
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